Our Workshop Target Region

To keep this workshop material concise and make it so you don’t need access to a server of cluster throughout the following demonstrations of KrigR, we will specify a set of locations in which we are interested.

The locations we focus on for this workshop are situated throughout eastern Germany and the north-western parts of the Czech Republic. Why do we focus on this particular part of the Earth? There are three reasons:

1. Topographical Heterogeneity - the area we select here contains large swaths of flat lowlands as well as some mountain ridges. This will make for visually pleasing plots and highlight the capability of kriging.
   
2. Geographic Scale - the area we are selecting here hits a certain sweet-spot for our purposes as its size makes it so that all KrigR functions run to completion in a relatively short time.
   
3. Familiarity - I was born and grew up in this region and have fond memories of the place. Please excuse my indulging in a bit of nostalgia.

Change the locations of interest at your own risk.

Using a different set of locations than the ones we specify here will change computational load and time as well as disk space required when working through the workshop material.

KrigR will be able to get you the data you want for the locations you desire, but computational requirements will vary.

Spatial Preferences in KrigR

KrigR is capable of learning about your spatial preferences in three ways:
1. As an extent input (a rectangular box).
2. As a SpatialPolygons input (a polygon or set of polygons).
3. As a set of locations stored in a data.frame.

To demonstrate the range of specifications permitted in KrigR, we make use of all three specifications. Masking out unnecessary areas from our analyses speeds up Kriging tremendously hence why we strongly suggest you make use of SpatialPolygons or data.frames whenever possible.

Extent_ext <- extent(c(9.87, 15.03, 49.89, 53.06))

Shape_shp <- ne_states(country = c("Germany", "Czech Republic"))
Shape_shp <- Shape_shp[Shape_shp$name_en %in% c("Saxony", "Saxony-Anhalt", "Thuringia", 
                                                "Ústí nad Labem Region", "Karlovy Vary Region"), ]

download.file("https://www.naturalearthdata.com/http//www.naturalearthdata.com/download/10m/cultural/ne_10m_admin_1_states_provinces.zip",
              destfile = "highres.zip")
unzip("highres.zip")
Shape_shp <- readOGR("ne_10m_admin_1_states_provinces.shp")
Shape_shp <- Shape_shp[Shape_shp$name_en %in% c("Saxony", "Saxony-Anhalt", "Thuringia",
                                                "Ústí nad Labem", "Karlovy Vary"), ]

load(file.path(Dir.Data, "Mountains_df.RData")) # load an sp object called Mountains_sp
Mountains_df

##          Mountain      Lon      Lat
## 1     Fichtelberg 12.95472 50.42861
## 2         Brocken 10.61722 51.80056
## 3 Großer Beerberg 10.74611 50.65944
## 4        Meluzína 13.00778 50.39028
## 5       Milešovka 13.93153 50.55523

Visualising our Study Setting

To finish our preparations for this workshop, let’s visualise the different locations of interest:

## Establish rectangular bounding box from extent
bbox <- as.numeric(as(Extent_ext, 'SpatialPolygons')@bbox)
names(bbox) <- c("left", "bottom", "right", "top")

## Make locations of mountains into SpatialPoints
Mountains_sp <- Mountains_df
coordinates(Mountains_sp) <- ~Lon+Lat

## download a map of the area specified by the extent
back_gg <- get_map(bbox, maptype = 'terrain')

## combine locations of interest into one plot
ggmap(back_gg, extent = "device") + # plot the extent area
  ## display the SpatialPolygons area
  geom_polygon(aes(x = long, y = lat, group = id), data = fortify(Shape_shp),
               colour = 'black', size = 1, fill = 'black', alpha = .5) + 
  ## add the data.frame data
  geom_point(aes(x = Lon, y = Lat), data = data.frame(Mountains_sp), 
             colour = "red", size = 4, pch = 13) + 
  ## some style additions
  theme_bw() + labs(x= "Longitude [°]", y = "Latitude  [°]") + 
  theme(plot.margin=unit(c(0, 1, 0, 1),"lines"))

In the above figure, the map area designates the extent specifications while the grey overlay display the SpatialPolygons preference and points of interest (form our data.frame input) are highlighted with red plotting symbols.

We are now ready to start the KrigR portion of the workshop!

Spatial Limitation

Let’s start with spatial limitation. As discussed previously, download_ERA() can handle a variety of inputs describing spatial preferences.

KrigR is capable of learning about your spatial preferences in three ways:
1. As an extent input (a rectangular box).
2. As a SpatialPolygons input (a polygon or set of polygons).
3. As a set of locations stored in a data.frame.
These spatial preferences are registered in KrigR functions using the Extent argument.

You might now ask yourself: How does KrigR achieve spatial limitation of the data? Couldn’t we just simply download only the data we are interested in?

The ECMWF CDS gives us tremendous capability of retrieving only the data we want. However, the CDS only recognises rectangular boxes (i.e., extents) for spatial limitation. Consequently, we always have to download data corresponding to a rectangular box in space. When informing KrigR of your spatial preferences using a data.frame or SpatialPolygons, download_ERA() automatically (1) identifies the smallest extent required by your input, (2) downloads data corresponding to this extent, and (3) masks our any data not queried by you.

Using KrigR’s spatial limitation features ensures faster computation and smaller file sizes (depending on file type).

In the following, I demonstrate how to use the Extent argument in download_ERA().

SpatialPolygonsRaw <- download_ERA(
  Variable = "2m_temperature",
  DataSet = "era5-land",
  DateStart = "1995-01-03",
  DateStop = "1995-01-03",
  TResolution = "day",
  TStep = 1,
  Extent = Shape_shp,
  Dir = Dir.Data,
  FileName = "SpatialPolygonsRaw",
  API_User = API_User,
  API_Key = API_Key
)
Plot_Raw(SpatialPolygonsRaw, 
         Shp = Shape_shp,
         Dates = c("01-1995", "02-1995", "03-1995", "04-1995")) 

Points_Raw <- download_ERA(
  Variable = "2m_temperature",
  DataSet = "era5-land",
  DateStart = "1995-01-01",
  DateStop = "1995-01-4",
  TResolution = "day",
  TStep = 1,
  Extent = Mountains_df,
  Buffer = 0.5,
  ID = "Mountain",
  Dir = Dir.Data,
  FileName = "PointsRaw",
  API_User = API_User,
  API_Key = API_Key
)
Plot_Raw(Points_Raw, Dates = c("01-1995", "02-1995", 
                              "03-1995", "04-1995")) + 
  geom_point(aes(x = Lon, y = Lat), data = data.frame(Mountains_sp), 
             colour = "green", size = 10, pch = 14)

Dynamical Data Uncertainty

With climate reanalyses, you also gain access to uncertainty flags of the data stored in the reanalysis product. For the ERA5-family of products, this uncertainty can be obtained by assessing the standard deviation of the 10 ensemble members which make up the underlying ERA5 model exercise.

For an aggregate understanding of data uncertainty, we also obtain dynamical uncertainty for our target region and time frame. For simplicity, we do so only for the SpatialPolygons specification.

With download_ERA() you can obtain this information as follows:

SpatialPolygonsEns <- download_ERA(
    Variable = "2m_temperature",
    DataSet = "era5",
    Type = "ensemble_members",
    DateStart = "1995-01-01",
    DateStop = "1995-01-04",
    TResolution = "day",
    TStep = 1,
    FUN = sd,
    Extent = Shape_shp,
    Dir = Dir.Data,
    FileName = "SpatialPolygonsEns",
    API_User = API_User,
    API_Key = API_Key
  )

Plot_Raw(SpatialPolygonsEns, Dates = c("01-1995", "02-1995", 
                                          "03-1995", "04-1995"),
         Shp = Shape_shp, COL = rev(viridis(100)))

As you can see here, there is substantial disagreement between the ensemble members of daily average temperatures across our study region. This uncertainty among ensemble members is greatest at high temporal resolution and becomes negligible at coarse temporal resolution. We document this phenomenon in this publication (Figure 1).

We will see how these uncertainties stack up against other sources of uncertainty when we arrive at aggregate uncertainty of our final product.

Covariates

First, we use the download_DEM() function which comes with KrigR to obtain elevation data as our covariate of choice. This produces two rasters:
1. A raster of training resolution which matches the input data in all attributes except for the data in each cell.
2. A raster of target resolution which matches the input data as closely as possible in all attributes except for the resolution (which is specified by the user).

Both of these products are bundled into a list where the first element corresponds to the training resolution and the second element contains the target resolution covariate data. Here, we specify a target resolution of .02.

This is how we specify download_DEM() to prepare DEM covariates for us:

Covs_ls <- download_DEM(Train_ras = SpatialPolygonsRaw, # the data we want to downscale
                        Target_res = .02, # the resolution we want to downscale to
                        Shape = Shape_shp, # extra spatial preferences
                        Dir = Dir.Covariates # where to store the covariate files
                        )

For now, let’s simply inspect our list of covariate rasters:

Covs_ls

## [[1]]
## class      : RasterLayer 
## dimensions : 34, 54, 1836  (nrow, ncol, ncell)
## resolution : 0.1000189, 0.09999998  (x, y)
## extent     : 9.726991, 15.12801, 49.75, 53.15  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : DEM 
## values     : 20.11554, 861.7248  (min, max)
## 
## 
## [[2]]
## class      : RasterLayer 
## dimensions : 204, 324, 66096  (nrow, ncol, ncell)
## resolution : 0.01666667, 0.01666667  (x, y)
## extent     : 9.72486, 15.12486, 49.74986, 53.14986  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : DEM 
## values     : 15.75, 1173  (min, max)

You will find that the target resolution covariate data comes at a resolution of 0.017 instead of the 0.02 resolution we specified. This happens because download_DEM() calls upon the raster::aggregate() function when aggregating the high-resolution covariate data to your desired target resolution and is thus only capable of creating target-resolution covariates in multiples of the base resolution of the GMTED 2010 DEM we are using as our default covariate. This happens only when the Target_res argument is specified to be a number.

Specifying the Target_res argument as a number will lead to best approximation of the desired resolution due to usage of the raster::aggregate() within download_DEM(). If you need an exact resolution to match pre-existing data, please refer to the third-party section of this workshop material.

Before moving on, let’s visualise the covariate data:

Plot_Covs(Covs_ls, Shape_shp)

Kriging

Kriging can be a very computationally expensive exercise.

The expense of kriging is largely determined by three factors:
1. Change in spatial resolution.
2. Number of cells containing data; i.e. Spatial Limitation.
3. Localisation of Kriging; i.e. Localisation of Results.

We explore two of these further down in this workshop material. For more information, please consult this publication (Figure 4).

Finally, we are ready to interpolate our input data given our covariates with the krigR() function:

SpatialPolygonsKrig <- krigR(Data = SpatialPolygonsRaw, # data we want to krig as a raster object
      Covariates_coarse = Covs_ls[[1]], # training covariate as a raster object
      Covariates_fine = Covs_ls[[2]], # target covariate as a raster object
      Keep_Temporary = FALSE, # we don't want to retain the individually kriged layers on our hard-drive
      Cores = 1, # we want to krig on just one core
      FileName = "SpatialPolygonsKrig", # the file name for our full kriging output
      Dir = Dir.Exports # which directory to save our final input in
      )

## Commencing Kriging

## Kriging of remaining 3 data layers should finish around: 2022-06-07 12:40:25

## 
  |                                                                                      
  |                                                                                |   0%
  |                                                                                      
  |====================                                                            |  25%
  |                                                                                      
  |========================================                                        |  50%
  |                                                                                      
  |============================================================                    |  75%
  |                                                                                      
  |================================================================================| 100%

Just like with the download_ERA() function, krigR() updates you on what it is currently working on. Again, I implemented this to make sure people don’t get too anxious staring at an empty console in R. If this feature is not appealing to you, you can turn this progress tracking off by setting verbose = FALSE in the function call to krigR().

For the rest of this workshop, I suppress messages from krigR() via other means so that when you execute, you get progress tracking.

There we go. As output of the krigR() function, we obtain a list of downscaled data as the first element and downscaling standard errors as the second list element. Let’s look at that:

SpatialPolygonsKrig[-3] # we talk later about the third element separately

## $Kriging_Output
## class      : RasterBrick 
## dimensions : 191, 309, 59019, 4  (nrow, ncol, ncell, nlayers)
## resolution : 0.01666667, 0.01666667  (x, y)
## extent     : 9.87486, 15.02486, 49.88319, 53.06653  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : var1.pred.1, var1.pred.2, var1.pred.3, var1.pred.4 
## min values :    269.5768,    266.5758,    265.7707,    261.4294 
## max values :    275.0120,    273.3299,    272.1343,    270.0602 
## 
## 
## $Kriging_SE
## class      : RasterBrick 
## dimensions : 191, 309, 59019, 4  (nrow, ncol, ncell, nlayers)
## resolution : 0.01666667, 0.01666667  (x, y)
## extent     : 9.87486, 15.02486, 49.88319, 53.06653  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : var1.stdev.1, var1.stdev.2, var1.stdev.3, var1.stdev.4 
## min values :    0.1251704,    0.1585938,    0.1427118,    0.1373274 
## max values :    0.1435884,    0.1745415,    0.1761995,    0.2816519

All the data has been downscaled and we do have uncertainties recorded for all of our outputs. Let’s visualise the data:

Plot_Krigs(SpatialPolygonsKrig, 
           Shp = Shape_shp,
           Dates = c("01-1995", "02-1995", "03-1995", "04-1995")
           )

As you can see, the elevation patterns show up clearly in our kriged air temperature output. Furthermore, you can see that our certainty of Kriging predictions drops on the 04/01/1995 in comparison to the preceding days. However, do keep in mind that a maximum standard error of 0.144, 0.175, 0.176, 0.282 (for each layer of our output respectively) on a total range of data of 5.435, 6.754, 6.364, 8.631 (again, for each layer in the output respectively) is evident of a downscaling result we can be confident in. We also demonstrated reliability of kriging in this publication (Figure 3).

Finally, this SpatialPolygons-informed downscaling took roughly 4 minutes on my machine (this may vary drastically on other devices).

Aggregate Uncertainty

Every climate data product is subject to an error-rate / range of data uncertainty. Unfortunately, almost none of the established climate data products communicate associated uncertainties. This leads to a dangerous overestimation of data credibility.

With the KrigR workflow, it is trivial to obtain uncertainty flags for all of your data - no matter the spatial or temporal resolution.

To understand the full certainty of our data obtained via the KrigR workflow, we should combine dynamical uncertainty with the statistical uncertainty we obtained from the krigR() function call above.

To do so, we require two data sets:
- SpatialPoylgonsKrig - created above containing statistical uncertainty in the second list position
- SpatialPoylgonsEns - created here; containing dynamical uncertainty

First, we assign them to objects with shorter names:

DynUnc <- SpatialPolygonsEns
KrigUnc <- SpatialPolygonsKrig[[2]]

Next, we need to align the rasters of statistical uncertainty (resolution: 0.017) and dynamical uncertainty (resolution: 0.5). As you can see, these are of differing resolutions and so cannot easily be combined using raster math. Instead, we first disaggregate the coarser-resolution raster (DynUnc) as disaggregation does not attempt any interpolation thus preserving the data, but representing it with smaller cells. To fix final remaining alignment issues, we allow for some resampling between the two raster:

EnsDisagg <- disaggregate(DynUnc, fact=res(DynUnc)[1]/res(KrigUnc)[1])
DynUnc <- resample(EnsDisagg, KrigUnc)

Finally, we combine the two uncertainty data products to form an aggregate uncertainty product:

FullUnc <- DynUnc + KrigUnc

Now, we are ready to plot our aggregate uncertainty:

Plot_Raw(FullUnc, 
         Shp = Shape_shp, 
         Dates = c("01-1995", "02-1995", "03-1995", "04-1995"), 
         COL = rev(viridis(100)))

As you can see, at short time-scales dynamic uncertainty eclipses statistical uncertainty. However, this phenomenon reverses at longer time-scales as shown in this publication (Figure 1).

Third-Party Data

Shape_shp <- ne_countries(country = "Germany")

## downloading gimms data
gimms_files <- downloadGimms(x = as.Date("2015-01-01"), # download from January 1982
                             y = as.Date("2015-12-31"), # download to December 1982
                             dsn = Dir.Data, # save downloads in data folder
                             quiet = FALSE # show download progress
                             )
## processing gimms data
gimms_raster <- rasterizeGimms(x = gimms_files, # the data we rasterize
                               remove_header = TRUE # we don't need the header of the data
                               )
indices <- monthlyIndices(gimms_files) # generate month indices from the data
gimms_raster_mvc <- monthlyComposite(gimms_raster, # the data
                                     indices = indices # the indices
                                     )
Negatives <- which(values(gimms_raster_mvc) < 0) # identify all negative values
values(gimms_raster_mvc)[Negatives] <- 0 # set threshold for barren land (NDVI<0)
gimms_raster_mvc <- crop(gimms_raster_mvc, extent(Shape_shp)) # crop to extent
gimms_mask <- KrigR::mask_Shape(gimms_raster_mvc[[1]], Shape = Shape_shp) # create mask ith KrigR-internal function to ensure all edge cells are contained
NDVI_ras <- mask(gimms_raster_mvc, gimms_mask) # mask out shape
NDVI_ras <- calc(NDVI_ras, fun = mean, na.rm = TRUE) # annual mean
writeRaster(NDVI_ras, format = "CDF", file = file.path(Dir.Data, "NDVI")) # save file

So what does this raster look like?

NDVI_ras

## class      : RasterStack 
## dimensions : 92, 108, 9936, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : 6, 15, 47.33333, 55  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :     layer 
## min values : 0.2430333 
## max values : 0.8339083

And a visualisation of the same:

Plot_Raw(NDVI_ras, 
         Shp = Shape_shp,
         Dates = "Mean NDVI 2015", 
         COL = viridis(100, begin = 0.5, direction = -1))

As stated above, we want to match this with our output.

KrigR Workflow

We could do this whole analysis in our three steps as outlined above, but why bother when the pipeline gets the job done just as well? The pipeline is a way by which you can specify download and downscaling of ERA5(-Land) data with GMTED 2010 DEM covariate data in just one function call.

Matching Kriging outputs with a pre-existing data set is as easy as plugging the pre-existing raster into the Target_res argument of the krigR() or the download_DEM() function.

This time we want to downscale from ERA5 resolution (roughly 30km) because the ERA5-Land data already matches the NDVI resolution (roughly 9km). Here’s how we do this:

NDVI_Krig <- krigR(
  ## download_ERA block
  Variable = '2m_temperature',
  Type = 'reanalysis',
  DataSet = 'era5',
  DateStart = '2015-01-01',
  DateStop = '2015-12-31',
  TResolution = 'year',
  TStep = 1,
  Extent = Shape_shp,
  API_User = API_User,
  API_Key = API_Key,
  SingularDL = TRUE,
  ## download_DEM block
  Target_res = NDVI_ras,
  Source = "Drive",
  ## krigR block
  Cores = 1,
  FileName = "AirTemp_NDVI.nc",
  nmax = 80, 
  Dir = Dir.Exports)

## download_ERA() is starting. Depending on your specifications, this can take a significant time.

## User 39340 for cds service added successfully in keychain

## Staging 1 download(s).

## Staging your request as a singular download now. This can take a long time due to size of required product.

## 0001_2m_temperature_2015-01-01_2015-12-31_year.nc download queried

## Requesting data to the cds service with username 39340

## - staging data transfer at url endpoint or request id:

##   a99a448b-7075-4bf7-9b60-bd52fb0731c1

## - timeout set to 10.0 hours

## - polling server for a data transfer \ polling server for a data transfer | polling
## server for a data transfer / polling server for a data transfer - polling server for a
## data transfer \ polling server for a data transfer | polling server for a data transfer /
## polling server for a data transfer - polling server for a data transfer \ polling server
## for a data transfer | polling server for a data transfer / polling server for a data
## transfer - polling server for a data transfer \ polling server for a data transfer |
## polling server for a data transfer / polling server for a data transfer - polling server
## for a data transfer \ polling server for a data transfer | polling server for a data
## transfer / polling server for a data transfer - polling server for a data transfer \
## polling server for a data transfer | polling server for a data transfer / polling server
## for a data transfer - polling server for a data transfer

## 
  |                                                                                      
  |                                                                                |   0%
  |                                                                                      
  |=============================================================                   |  76%
  |                                                                                      
  |================================================================================| 100%

## - moved temporary file to -> C:/Users/erike/Documents/Work/Z - Conferences-Workshops_Presentations/2022_06_15-17 - [NordOikos]/KrigR Workshop/Exports/0001_2m_temperature_2015-01-01_2015-12-31_year.nc
## - request purged from queue!
## Checking for known data issues.
## Loading downloaded data for masking and aggregation.
## Masking according to shape/buffer polygon
## Aggregating to temporal resolution of choice

## 
  |                                                                                      
  |                                                                                |   0%
  |                                                                                      
  |===========================                                                     |  33%
  |                                                                                      
  |=====================================================                           |  67%
  |                                                                                      
  |================================================================================| 100%
## 

## Commencing Kriging

## 
  |                                                                                      
  |                                                                                |   0%
  |                                                                                      
  |================================================================================| 100%

As you can see, executing krigR() as a pipeline updates you on every single step along the way.

So? Did we match the pre-existing data?

NDVI_Krig[[1]]

## class      : RasterBrick 
## dimensions : 92, 108, 9936, 1  (nrow, ncol, ncell, nlayers)
## resolution : 0.08333333, 0.08333333  (x, y)
## extent     : 6, 15, 47.33333, 55  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## source     : memory
## names      : var1.pred 
## min values :  275.9705 
## max values :  285.7357

We nailed this!

Let’s take one final look at our (A) raw ERA5 data, (B) NDVI data, (C) Kriged ERA5 data, and (D) standard error of our Kriging output:

### ERA-Plot
ERA_df <- as.data.frame(raster(file.path(Dir.Exports, "2m_temperature_2015-01-01_2015-12-31_year.nc")), xy = TRUE) # turn raster into dataframe
colnames(ERA_df)[c(-1,-2)] <- "Air Temperature 2015 (ERA5)"
ERA_df <- gather(data = ERA_df, key = Values, value = "value", colnames(ERA_df)[c(-1,-2)]) #  make ggplot-ready
Raw_plot <- ggplot() + # create a plot
  geom_raster(data = ERA_df , aes(x = x, y = y, fill = value)) + # plot the raw data
  facet_wrap(~Values) + # split raster layers up
  theme_bw() + labs(x = "Longitude", y = "Latitude") + # make plot more readable
  scale_fill_gradientn(name = "Air Temperature [K]", colours = inferno(100)) + # add colour and legend
  geom_polygon(data = Shape_shp, aes(x = long, y = lat, group = group), colour = "black", fill = "NA") # add shape
### NDVI-Plot
NDVI_df <- as.data.frame(NDVI_ras, xy = TRUE) # turn raster into dataframe
colnames(NDVI_df)[c(-1,-2)] <- "NDVI 2015"
NDVI_df <- gather(data = NDVI_df, key = Values, value = "value", colnames(NDVI_df)[c(-1,-2)]) #  make ggplot-ready
NDVI_plot <- ggplot() + # create a plot
  geom_raster(data = NDVI_df , aes(x = x, y = y, fill = value)) + # plot the raw data
  facet_wrap(~Values) + # split raster layers up
  theme_bw() + labs(x = "Longitude", y = "Latitude") + # make plot more readable
  scale_fill_gradientn(name = "NDVI", colours = rev(terrain.colors(100))) + # add colour and legend
  geom_polygon(data = Shape_shp, aes(x = long, y = lat, group = group), colour = "black", fill = "NA") # add shape
### KRIGED-Plots
Dates = c("Kriged Air Temperature 2015 (NDVI Resolution)")
Type_vec <- c("Prediction", "Standard Error") # these are the output types of krigR
Colours_ls <- list(inferno(100), rev(viridis(100))) # we want separate colours for the types
Plots_ls <- as.list(NA, NA) # this list will be filled with the output plots
KrigDF_ls <- as.list(NA, NA) # this list will be filled with the output data
for(Plot in 1:2){ # loop over both output types
  Krig_df <- as.data.frame(NDVI_Krig[[Plot]], xy = TRUE) # turn raster into dataframe
  colnames(Krig_df)[c(-1,-2)] <- paste(Type_vec[Plot], Dates) # set colnames
  Krig_df <- gather(data = Krig_df, key = Values, value = "value", colnames(Krig_df)[c(-1,-2)]) # make ggplot-ready
  Plots_ls[[Plot]] <- ggplot() + # create plot
    geom_raster(data = Krig_df , aes(x = x, y = y, fill = value)) + # plot the kriged data
    facet_wrap(~Values) + # split raster layers up
    theme_bw() + labs(x = "Longitude", y = "Latitude") + # make plot more readable
    scale_fill_gradientn(name = "Air Temperature [K]", colours = Colours_ls[[Plot]]) + # add colour and legend
    theme(plot.margin = unit(c(0, 0, 0, 0), "cm")) + # reduce margins (for fusing of plots)
    geom_polygon(data = Shape_shp, aes(x = long, y = lat, group = group), colour = "black", fill = "NA") # add shape
  KrigDF_ls[[Plot]] <- Krig_df
} # end of type-loop

plot_grid(plotlist = list(Raw_plot, NDVI_plot, Plots_ls[[1]], Plots_ls[[2]]), 
          nrow = 2, labels = "AUTO")

So what can we learn from this? Let’s plot the relation between temperature and NDVI:

plot_df <- as.data.frame(cbind(KrigDF_ls[[1]][,4], 
                               KrigDF_ls[[2]][,4],
                               NDVI_df[,4]))
colnames(plot_df) <- c("Temperature", "Uncertainty", "NDVI")
ggplot(plot_df,
       aes(x = Temperature, y = NDVI, size = Uncertainty)) + 
  geom_point(alpha = 0.15) + 
  theme_bw()

Looks like NDVI increases as mean annual temperatures rise, but reaches a peak around 281-282 Kelvin with a subsequent decrease as mean annual temperatures rise higher.

Third-Party Data Covariates

In this publication, we demonstrate how soil moisture data can be reliably statistically downscaled using soil property data which we obtain from the Land-Atmosphere Interaction Research Group at Sun Yat-sen University.

Below, you will find the code needed to obtain the data of global coverage at roughly 1km spatial resolution. The code chunk below also crops and masks the data according to our study region and subsequently deletes the storage-heavy global files (3.5GB each in size). This process takes a long time due to download speeds.

# documentation of these can be found here http://globalchange.bnu.edu.cn/research/soil4.jsp
SoilCovs_vec <- c("tkdry", "tksat", "csol", "k_s", "lambda", "psi", "theta_s") # need these names for addressing soil covariates
if(!file.exists(file.path(Dir.Covariates, "SoilCovs.nc"))){
  print("#### Loading SOIL PROPERTY covariate data. ####") 
  # create lists to combine soil data into one
  SoilCovs_ls <- as.list(rep(NA, length(SoilCovs_vec)))
  names(SoilCovs_ls) <- c(SoilCovs_vec)
  ## Downloading, unpacking, and loading
  for(Soil_Iter in SoilCovs_vec){
    if(!file.exists(file.path(Dir.Covariates, paste0(Soil_Iter, ".nc")))) { # if not downloaded and processed yet
      print(paste("Handling", Soil_Iter, "data."))
      Dir.Soil <- file.path(Dir.Covariates, Soil_Iter)
      dir.create(Dir.Soil)
      download.file(paste0("http://globalchange.bnu.edu.cn/download/data/worldptf/", Soil_Iter,".zip"),
                    destfile = file.path(Dir.Soil, paste0(Soil_Iter, ".zip"))
      ) # download data
      unzip(file.path(Dir.Soil, paste0(Soil_Iter, ".zip")), exdir = Dir.Soil) # unzip data
      File <- list.files(Dir.Soil, pattern = ".nc")[1] # only keep first soil layer
      Soil_ras <- raster(file.path(Dir.Soil, File)) # load data
      SoilCovs_ls[[which(names(SoilCovs_ls) == Soil_Iter)]] <- Soil_ras # save to list
      writeRaster(x = Soil_ras, filename = file.path(Dir.Covariates, Soil_Iter), format = "CDF")
      unlink(Dir.Soil, recursive = TRUE)
    }else{
      print(paste(Soil_Iter, "already downloaded and processed."))
      SoilCovs_ls[[which(names(SoilCovs_ls) == Soil_Iter)]] <- raster(file.path(Dir.Covariates, paste0(Soil_Iter, ".nc")))
    }
  }
  ## data handling and manipulation
  SoilCovs_stack <- stack(SoilCovs_ls) # stacking raster layers from list
  SoilCovs_stack <- crop(SoilCovs_stack, extent(BCq_ras)) # cropping to extent of data we have
  SoilCovs_mask <- KrigR::mask_Shape(SoilCovs_stack[[1]], Shape = Shape_shp) # create mask with KrigR-internal function to ensure all edge cells are contained
  SoilCovs_stack <- mask(SoilCovs_stack, SoilCovs_mask) # mask out shape
  ## writing the data
  writeRaster(x = SoilCovs_stack, filename = file.path(Dir.Covariates, "SoilCovs"), format = "CDF")
  ## removing the global files due to their size
  unlink(file.path(Dir.Covariates, paste0(SoilCovs_vec, ".nc")))
}

Let’s have a look at these data:

SoilCovs_stack <- stack(file.path(Dir.Covariates, "SoilCovs.nc"))
names(SoilCovs_stack) <- SoilCovs_vec
SoilCovs_stack

## class      : RasterStack 
## dimensions : 408, 648, 264384, 7  (nrow, ncol, ncell, nlayers)
## resolution : 0.008333333, 0.008333333  (x, y)
## extent     : 9.725, 15.125, 49.75, 53.15  (xmin, xmax, ymin, ymax)
## crs        : +proj=longlat +datum=WGS84 +no_defs 
## names      :         tkdry,         tksat,          csol,           k_s,        lambda,           psi,       theta_s 
## min values :  5.200000e-02,  1.337000e+00,  2.141000e+06,  5.212523e+00,  8.600000e-02, -5.307258e+01,  3.230000e-01 
## max values :  2.070000e-01,  2.862000e+00,  2.346400e+06,  2.461686e+02,  3.330000e-01, -5.205317e+00,  5.320000e-01

Now we need to establish target and training resolution of our covariate data.

First, we focus on the training resolution covariate data. We match our covariate data to our spatial product which we wish to downscale by resampling the covariate data to the coarser resolution:

Coarsecovs <- resample(x = SoilCovs_stack, y = BCq_ras)

Second, we aggregate the covariate data to our desired resolution. In this case, 0.02 as done previously:

Finecovs <- aggregate(SoilCovs_stack, fact = 0.02/res(SoilCovs_stack)[1])

Finally, we combine these into a list like the output of download_DEM():

Covs_ls <- list(Coarsecovs, Finecovs)
Plot_Covs(Covs = Covs_ls, Shape_shp)

Our development goals include creating a function that automatically carries out all of the above for you with a specification alike to download_DEM().

Kriging Third-Party Data

Finally, we can statistically downscale our soil moisture data using the soil property covariates. For this, we need to specify a new KrigingEquation.

With the KrigingEquation argument, you may specify non-linear combinations of covariates for your call to krigR().

If you don’t specify a KrigingEquation in krigR() and your covariates do not contain a layer called "DEM", krigR() will notify you that its default formula cannot be executed and will attempt to build an additive formula from the data it can find. krigr() will inform you of this and ask for your approval before proceeding.

This auto-generated formula would be the same as the one we specify here - an additive combination of all covariates found both at coarse and fine resolutions. Of course, this formula can also be specified to reflect interactive effects.

Here, I automate the generation of our KrigingEquation:

KrigingEquation <- paste0("ERA ~ ", paste(SoilCovs_vec, collapse = " + "))
KrigingEquation

## [1] "ERA ~ tkdry + tksat + csol + k_s + lambda + psi + theta_s"

I only hand the mean annual soil moisture layer to the kriging call and I set nmax to 80 to approximate a typical weather system in size:

BC_Water_Krig  <- krigR(Data = BCq_ras[[12]], 
                    Covariates_coarse = Covs_ls[[1]], 
                    Covariates_fine = Covs_ls[[2]],
                    KrigingEquation = KrigingEquation, 
                    FileName = "BC_Water_Krig",
                    Dir = Dir.Exports,
                    nmax = 80
)

Interpolating this data is just like statistically downscaling any other soil moisture product and can be done without any problems.

Look at how well the river Elbe sows up in this!

Plot_Krigs(BC_Water_Krig[-3], 
           Shp = Shape_shp,
           Dates = "BIO12 - Annual Mean Soil Moisture",
           columns = 2)

Obtaining ERA5(-Land) Data

Now, let’s obtain the historical baseline from ERA5-Land for the same time-period as our CMIP6 historical data.

if(!file.exists(file.path(Dir.Data, "Germany_Hist_ERA.nc"))){
  Hist_ERA_ras <- download_ERA(Variable = "2m_temperature",
                               DateStart = "1981-01-01",
                               DateStop = "1999-12-31",
                               TResolution = "month",
                               TStep = 1,
                               Extent = Shape_shp,
                               Dir = Dir.Data,
                               FileName = "Germany_Hist_ERA", 
                               API_Key = API_Key,
                               API_User = API_User,
                               SingularDL = TRUE)
  Index <- rep(1:12, length = nlayers(Hist_ERA_ras))
  Hist_ERA_ras <- stackApply(Hist_ERA_ras, indices = Index, fun = mean)
  writeRaster(Hist_ERA_ras, filename = file.path(Dir.Data, "Germany_Hist_ERA"), format = "CDF")
}
Hist_ERA_ras <- mean(stack(file.path(Dir.Data, "Germany_Hist_ERA.nc")))

Obtaining Projection Data

Here, we use CMIP6 projection data manually sourced from the ECMWF CDS distribution.

Our development goals include development of download_ERA() to work with other ECWMF CDS data sets aside from ERA5(-Land). This includes this CMIP6 data set.

train_HIST <- mean(stack(file.path(Dir.Data, "historical_tas_1981-2000.nc")))
train_HIST <- crop(train_HIST,extent(Hist_ERA_ras))
train_mask <- KrigR::mask_Shape(train_HIST, Shape_shp)
train_HIST <- mask(train_HIST, train_mask)

train_SSP <- mean(stack(file.path(Dir.Data, "ssp585_tas_2041-2060.nc")))
train_SSP <- crop(train_SSP,extent(Hist_ERA_ras))
train_mask <- KrigR::mask_Shape(train_SSP, Shape_shp)
train_SSP <- mask(train_SSP, train_mask)

Plot_Raw(stack(train_HIST, train_SSP), 
         Shp = Shape_shp,
         Dates = c("Historic CMIP6", "Future CMIP6"))

Establishing Kriged Products

For the first time in this workshop material, we will push our spatial resolution to the finest scale supported by our default GMTED 2010 DEM covariate data: 0.008333 / ~1km.

These operations take quite some time - grab a tea or coffee, go for a walk, or stretch a bit.

The downscaling calls should be familiar by now so I will forego explaining them. In case, the following code snippets do not make sense to you, please consult the portion of this workshop concerned with statistical downscaling.

## Covariate Data
GMTED_DE <- download_DEM(
  Train_ras = train_HIST,
  Target_res = 0.008334,
  Shape = Shape_shp,
  Keep_Temporary = TRUE,
  Dir = Dir.Covariates
)
## Kriging
Output_HIST <- krigR(
  Data = train_HIST,
  Covariates_coarse = GMTED_DE[[1]], 
  Covariates_fine = GMTED_DE[[2]],  
  Keep_Temporary = FALSE,
  Cores = 1,
  Dir = Dir.Exports,  
  FileName = "DE_CMIP-HIST", 
  nmax = 40
)

Plot_Krigs(Output_HIST,
           Shp = Shape_shp,
           Dates = "CMIP6 Historical", columns = 2)

## Covariate Data
GMTED_DE <- download_DEM(
  Train_ras = train_SSP,
  Target_res = 0.008334,
  Shape = Shape_shp,
  Keep_Temporary = TRUE,
  Dir = Dir.Covariates
)
## Kriging
Output_SSP <- krigR(
  Data = train_SSP,
  Covariates_coarse = GMTED_DE[[1]], 
  Covariates_fine = GMTED_DE[[2]],   
  Keep_Temporary = FALSE,
  Cores = 1,
  Dir = Dir.Exports,  
  FileName = "DE_SSP585_2041-2060", 
  nmax = 40
)

Plot_Krigs(Output_SSP,
           Shp = Shape_shp,
           Dates = "CMIP6 Future", columns = 2)

## Covariate Data
GMTED_DE <- download_DEM(
  Train_ras = Hist_ERA_ras,
  Target_res = 0.008334,
  Shape = Shape_shp,
  Keep_Temporary = TRUE,
  Dir = Dir.Covariates
)
## Kriging
Output_ERA <- krigR(
  Data = Hist_ERA_ras,
  Covariates_coarse = GMTED_DE[[1]], 
  Covariates_fine = GMTED_DE[[2]],   
  Keep_Temporary = FALSE,
  Cores = 1,
  Dir = Dir.Exports,  
  FileName = "DE_hist", 
  nmax = 40
)

Plot_Krigs(Output_ERA,
           Shp = Shape_shp,
           Dates = "ERA5-Land Historical", columns = 2)

Putting It All Together

To establish a final product of high-resolution climate projection data, we simply add the difference between the kriged CMIP6 products to the kriged ERA5-Land product:

## Creating Difference and Projection raster
Difference_ras <- Output_SSP[[1]] - Output_HIST[[1]]
Projection_ras <- Output_ERA[[1]] + Difference_ras
## Adding min and max values to ocean cells to ensure same colour scale
Output_ERA[[1]][10] <- maxValue(Projection_ras)
Output_ERA[[1]][12] <- minValue(Projection_ras)
Projection_ras[10] <- maxValue(Output_ERA[[1]])
Projection_ras[12] <- minValue(Output_ERA[[1]])
## Individual plots
A_gg <- Plot_Raw(Output_ERA[[1]], Shp = Shape_shp, 
                 Dates = "Historical ERA5-Land (1981-2000)")
B_gg <- Plot_Raw(Difference_ras[[1]], Shp = Shape_shp, 
                 Dates = "Anomalies of SSP585 - Historical CMIP-6",
                 COL = rev(viridis(100)))
C_gg <- Plot_Raw(Projection_ras[[1]], Shp = Shape_shp, 
                 Dates = "Future Projection (ERA5-Land + Anomalies)")
## Fuse the plots into one big plot
ggPlot <- plot_grid(plotlist = list(A_gg, B_gg, C_gg), 
                    ncol = 3, labels = "AUTO") 
ggPlot